NATRE region variance budget#

Reproducing the Ferrari & Polzin (2005) estimate using NATRE data, and comparing differences with our version

TODO#

  1. [ ] Microscale variance budget

    1. Could estimate \(⟨u_t T_t⟩ ⋅ ∇T_m = -⟨K_t T_z⟩ T_z^m\) term (monthly mean heat flux X monthly mean vertical gradient) using the model; just like I do with the microstructure data.

    2. Then add the mesoscale dissipation term to get total \(⟨χ⟩\)

Setup#

Hide code cell source 
%load_ext watermark
%load_ext autoreload
%matplotlib inline
%matplotlib inline

import cf_xarray as cfxr  # noqa
import cf_xarray.units  # noqa
import colorcet
import dcpy
import distributed
import gsw
import holoviews as hv
import hvplot.xarray  # noqa
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import pint_xarray  # noqa
import scipy as sp
import tqdm

hv.notebook_extension("bokeh")
import xgcm
from datatree import DataTree
from IPython.display import Image

import xarray as xr
from xarray.tests import raise_if_dask_computes

%aimport eddydiff
ed = eddydiff

xr.set_options(keep_attrs=True)

plt.rcParams["figure.dpi"] = 140
plt.rcParams["savefig.dpi"] = 200
plt.style.use("ggplot")
plt.figure()
plt.close()

%watermark -iv
xr.DataArray([1.0])

holoviews  : 1.15.4
dcpy       : 0.2.1.dev45+g67ec22a.d20230504
colorcet   : 3.0.1
pint_xarray: 0.3
tqdm       : 4.65.0
distributed: 2023.3.2
numpy      : 1.23.5
hvplot     : 0.8.3
matplotlib : 3.7.1
xarray     : 2023.3.0
eddydiff   : 0.1
gsw        : 3.6.16
xgcm       : 0.8.1
scipy      : 1.10.1
cf_xarray  : 0.8.1.dev0+gce9cd5a.d20230327

<xarray.DataArray (dim_0: 1)>
array([1.])
Dimensions without coordinates: dim_0
Hide code cell source 
if "client" in locals():
    client.cluster.close()
    client.close()
client = distributed.Client(n_workers=3, threads_per_worker=1)
client

Client

Client-868c6dde-150a-11ee-ba4b-acde48001122

Connection method: Cluster object Cluster type: distributed.LocalCluster
Dashboard: http://127.0.0.1:8787/status

Cluster Info

 
Image("../images/natre-large-scale.png")
_images/3404b5a3dd2cc1c8ee8b8c190f8e0940977ec605476088f48503cef1bc59deba.png

Prepare various production and dissipation terms#

Micro#

a05_binned = xr.load_dataset("../datasets/a05-2015.nc")

micro = xr.load_dataset("../datasets/natre-density-binned.nc")
Hide code cell source 
natre_fine = xr.load_dataset("../datasets/natre-finescale.nc")
natre_fine_avg = (
    natre_fine.groupby_bins("pressure", bins=np.arange(250, 2001, 200))
    .mean()
    .cf.guess_coord_axis()
)
natre_fine_avg.pressure_bins.attrs = {
    "positive": "down",
    "standard_name": "sea_water_pressure",
}

Argo#

Hide code cell source 
argo_fine_profiles = xr.load_dataset(
    "../datasets/argomix_natre.zarr", engine="zarr"
).drop("DIRECTION")
argo_fine_profiles["Kt"] = argo_fine_profiles["Kρ"]
argo_fine_profiles["Kt"].attrs["long_name"] = "$K_T$"
# argo_fine_profiles = argo_fine_profiles.where(argo_fine_profiles.ε < 0.1)
bins = cfxr.bounds_to_vertices(micro.gamma_n_bounds, "bound").data
argo_fine_dens = (
    argo_fine_profiles.rename({"γmean": "gamma_n"})
    .reset_coords("pressure")
    .groupby_bins("gamma_n", bins=bins)
    .mean()
)
argo_fine_dens["gamma_n_bins"].attrs = {
    "standard_name": "neutral_density",
    "axis": "Z",
    "positive": "down",
}
argo_fine_dens["npts"] = (
    argo_fine_profiles.χ.rename({"γmean": "gamma_n"})
    .groupby_bins("gamma_n", bins=bins)
    .count()
)
argo_fine_dens = argo_fine_dens.set_coords("pressure")

argo_fine_depth = argo_fine_profiles.groupby_bins(
    "pressure", np.arange(250, 2050, 200)
).mean()
argo_fine_depth["pressure_bins"].attrs = {"standard_name": "sea_water_pressure"}

ECCO#

ecco_tile2 = xr.open_zarr("../datasets/ecco-tile2-sigma2.zarr")
ecco_tile2["THETA"].attrs["coordinates"] = "XC YC z_σ"
ecco_tile2["SALT"].attrs["coordinates"] = "XC YC z_σ"

ecco_natre = ecco_tile2.isel(i=slice(10, 18), j=slice(10, 20)).load()
ecco_natre_profile = ed.pop.calc_mean_redivar_profile(ecco_natre.mean("time"))
ecco_natre_initial = ed.pop.calc_mean_redivar_profile(ecco_natre.isel(time=0))
ecco_natre_final = ed.pop.calc_mean_redivar_profile(ecco_natre.isel(time=-1))
plt.legend(["mean", "initial", "final"])
ecco_natre_profile

<xarray.Dataset>
Dimensions:      (σ: 60, bounds: 2)
Coordinates:
    tile         int64 2
    z_σ          (σ) float32 -64.37 -64.77 -65.33 ... -4.689e+03 -4.689e+03
  * σ            (σ) float32 34.15 34.15 34.17 34.18 ... 37.2 37.2 37.2 37.2
    z_σ_bounds   (σ, bounds) float32 -64.18 -64.57 ... -4.689e+03 -4.689e+03
Dimensions without coordinates: bounds
Data variables:
    SALT         (σ) float32 36.94 36.94 36.94 36.93 36.93 ... nan nan nan nan
    THETA        (σ) float32 21.72 21.69 21.65 21.6 21.47 ... nan nan nan nan
    delT2_plane  (σ) float64 1.35e-12 1.339e-12 1.326e-12 ... 0.0 0.0 0.0
sub = ecco_tile2.SALT.sel(σ=36, method="nearest").load()

sub.sel(i=slice(40), j=slice(40)).hvplot.quadmesh(
    geo=True,
    widget_type="scrubber",
    clim=(35.3, 35.8),
    cmap="bmy_r",
    widget_location="bottom",
)

ecco_natre_initial.delT2_plane.plot()

[<matplotlib.lines.Line2D at 0x17473b940>]

ecco = (
    xr.open_dataset("../datasets/ecco-chi.nc")
    .sel(time=slice("2000", None))
    .mean("time")
)
for var in ecco.cf.coordinates["vertical"]:
    ecco[var] = -1 * ecco[var]
kredi = ed.format_ecco2(
    xr.open_dataset("~/datasets/ecco/interp_Ks_ECCOv4.nc")
).Ks_ECCOv4.rename({"pres": "k"})
ecco["Ke"] = (
    kredi.sel(lat=slice(22, 28), lon=slice(322, 336)).mean(["lat", "lon"]).variable
)

ecco_natre_profile.update(ecco)
ecco_natre_profile["z_σ"] = -1 * ecco_natre_profile.z_σ
del ecco_natre_profile["σ"].attrs["positive"]
ecco_natre_profile["RediVar"] = ecco_natre_profile.GM_CHI
ecco_natre_profile

/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/xarray/core/indexing.py:661: RuntimeWarning: deallocating CachingFileManager(<class 'netCDF4._netCDF4.Dataset'>, '/Users/dcherian/work/eddydiff/datasets/ecco-chi.nc', mode='r', kwargs={'clobber': True, 'diskless': False, 'persist': False, 'format': 'NETCDF4'}, manager_id='85ce065b-6ab5-4525-ac28-ceecad6c8c97'), but file is not already closed. This may indicate a bug.
  self.array = NumpyIndexingAdapter(np.asarray(self.array))

<xarray.Dataset>
Dimensions:      (σ: 60, node: 2, k: 50, bounds: 2, k_u: 50, k_l: 50, k_p1: 51)
Coordinates: (12/20)
    tile         int64 2
    z_σ          (σ) float32 64.37 64.77 65.33 ... 4.689e+03 4.689e+03 4.689e+03
  * σ            (σ) float32 34.15 34.15 34.17 34.18 ... 37.2 37.2 37.2 37.2
    z_σ_bounds   (σ, bounds) float32 -64.18 -64.57 ... -4.689e+03 -4.689e+03
  * k            (k) int64 0 1 2 3 4 5 6 7 8 9 ... 40 41 42 43 44 45 46 47 48 49
  * k_u          (k_u) int64 0 1 2 3 4 5 6 7 8 9 ... 41 42 43 44 45 46 47 48 49
    ...           ...
    drF          (k) float32 10.0 10.0 10.0 10.0 ... 387.5 410.5 433.5 456.5
    PHrefC       (k) float32 49.05 147.1 245.2 ... 4.944e+04 5.357e+04 5.794e+04
    PHrefF       (k_p1) float32 0.0 98.1 196.2 ... 5.145e+04 5.57e+04 6.018e+04
    rhoRef       (k) float32 1.024e+03 1.024e+03 ... 1.052e+03 1.054e+03
    drC_kl       (k_l) float32 5.0 10.0 10.0 10.0 ... 376.0 399.0 422.0 445.0
  * node         (node) object 'adj' 'const'
Dimensions without coordinates: bounds
Data variables: (12/23)
    SALT         (σ) float32 36.94 36.94 36.94 36.93 36.93 ... nan nan nan nan
    THETA        (node, k) float32 22.61 22.59 22.54 22.38 ... 1.021 0.2023 0.0
    delT2_plane  (σ) float64 1.35e-12 1.339e-12 1.326e-12 ... 0.0 0.0 0.0
    GM_CHIXO     (node, k) float32 -1.491e-09 -2.558e-09 ... -7.767e-14 0.0
    GM_CHIZ      (node, k_l) float32 0.0 4.938e-09 5.708e-09 ... 1.499e-14 0.0
    GM_ubT       (node, k) float32 4.472e+05 -5.016e+04 ... 6.831e+03 0.0
    ...           ...
    Tx           (node, k) float32 -1.636e-06 -1.639e-06 ... -1.105e-06 0.0
    Ty           (node, k) float32 -1.341e-06 -1.341e-06 ... -1.417e-08 0.0
    Tz           (node, k_l) float32 nan 0.001619 ... 0.001941 0.0004545
    GM_CHI       (node, k) float32 1.014e-08 9.531e-09 ... 4.292e-08 nan
    Ke           (k) float32 1.85e+03 2.002e+03 2.06e+03 ... 984.2 958.3 948.3
    RediVar      (node, k) float32 1.014e-08 9.531e-09 ... 4.292e-08 nan

POP#

POP 1/10° Variance budget#

Hide code cell source 
pop_hires_natre = (
    xr.load_dataset("../datasets/natre-pop-variance.nc")
    .pint.quantify()
    .pint.to("K**2/s")
)
pop_hires_natre

<xarray.Dataset>
Dimensions:  (depth: 62)
Coordinates:
  * depth    (depth) float64 5.0 15.0 25.0 ... 5.375e+03 5.625e+03 5.875e+03
Data variables:
    BC       (depth) float64 [K²/s] 3.765e-09 2.352e-09 ... -3.835e-12
    HDIFF    (depth) float64 [K²/s] -3.8e-10 -3.846e-10 ... -3.298e-14
    PKC      (depth) float64 [K²/s] 4.486e-10 6.327e-11 ... 4.232e-15 0.0
    VMIX     (depth) float64 [K²/s] -9.317e-09 -1.392e-09 ... 4.307e-14
    DISS     (depth) float64 [K²/s] -9.697e-09 -1.777e-09 ... 1.01e-14
pop_hires_natre[["HDIFF", "VMIX"]].to_array().plot.line(
    hue="variable", y="depth", ylim=(3000, 0)
)

[<matplotlib.lines.Line2D at 0x178a4ddf0>,
 <matplotlib.lines.Line2D at 0x178a4de50>]

POP 1/10° state vars climatology#

pop_hires_clim = (
    xr.load_dataset("../datasets/pop-hires-natre-coarsened-annual-climatology.nc")
    .pint.quantify()
    .squeeze("cycle")
    .map(dcpy.util.to_base_units)
)
# pop_hires_clim['z_σ'] = pop_hires_clim.z_σ.cf.ffill("Z")
# pop_hires_clim = pop_hires_clim.cf.add_bounds("z_σ", dim="Z").isel(σ=slice(-7))
pop_hires_clim["TEMP"].attrs["standard_name"] = "sea_water_potential_temperature"
# NOTE: This only adds delT2, and averages in box
pop_hires_clim_profile = ed.pop.calc_mean_redivar_profile(pop_hires_clim)

# interpolate to pop_natre Z grid so we can estimate Ke
pop_hires_clim_profile = (
    pop_hires_clim_profile.swap_dims({"σ": "z_σ"})
    .drop_vars("z_σ_bounds")
    .dropna("z_σ")
    .interp(z_σ=pop_hires_natre.depth)
    .cf.add_bounds("depth")
)
pop_hires_clim_profile["Ke"] = (
    -1 * pop_hires_natre.DISS
) / pop_hires_clim_profile.delT2_plane.where(pop_hires_clim_profile.delT2_plane > 5e-13)

pop_hires_clim_profile

/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/xarray/core/duck_array_ops.py:675: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
  return push(array, n, axis)
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/scipy/interpolate/_polyint.py:120: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
  yi = np.asarray(yi)
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/scipy/interpolate/_interpolate.py:508: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
  y = array(y, copy=self.copy)

<xarray.Dataset>
Dimensions:       (depth: 62, bounds: 2)
Coordinates:
    time          object 0052-01-01 05:17:48.750000
    σ             (depth) float64 nan nan nan nan nan ... 37.2 nan nan nan nan
    z_σ           (depth) float64 5.0 15.0 25.0 ... 5.625e+03 5.875e+03
  * depth         (depth) float64 5.0 15.0 25.0 ... 5.625e+03 5.875e+03
    depth_bounds  (depth, bounds) float64 0.0 10.0 10.0 ... 5.75e+03 6e+03
    z_σ_bounds    (depth, bounds) float64 0.0 10.0 10.0 ... 5.75e+03 6e+03
Dimensions without coordinates: bounds
Data variables:
    TEMP          (depth) float64 nan nan nan nan nan ... 275.0 nan nan nan nan
    SALT          (depth) float64 nan nan nan nan nan ... nan nan nan nan
    delT2         (depth) float64 nan nan nan nan nan ... nan nan nan nan
    delT2_plane   (depth) float64 nan nan nan nan nan ... nan nan nan nan
    Ke            (depth) float64 [K²/s] nan nan nan nan nan ... nan nan nan nan
Attributes: (12/14)
    history:                    Mon Apr 15 11:06:38 2019: ncks -A 1st_half/ta...
    title:                      g.e20.G.TL319_t13.control.001_hfreq
    model_doi_url:              https://doi.org/10.5065/D67H1H0V
    time_period_freq:           day_5
    Conventions:                CF-1.0; http://www.cgd.ucar.edu/cms/eaton/net...
    yrs_averaged:               42-61
    ...                         ...
    calendar:                   All years have exactly  365 days.
    start_time:                 This dataset was created on 2019-01-16 at 20:...
    contents:                   Diagnostic and Prognostic Variables
    revision:                   $Id: tavg.F90 89091 2018-04-30 15:58:32Z altu...
    history_of_appended_files:  Mon Apr 15 11:06:38 2019: Appended file 1st_h...
    NCO:                        netCDF Operators version 4.7.4 (http://nco.sf...

POP 1°#

pop_1deg_natre_profile = xr.load_dataset(
    "../datasets/pop-1deg-redivar-natre.nc"
).pint.quantify()
pop_1deg_natre_profile.σ.attrs["axis"] = "Z"
pop_1deg_natre_profile["z_σ"] = pop_1deg_natre_profile.z_σ.cf.ffill("Z")
pop_1deg_natre_profile = pop_1deg_natre_profile.cf.add_bounds("z_σ", dim="Z")
pop_1deg_natre_profile["delT2_plane"] = pop_1deg_natre_profile.delT2
pop_1deg_natre_profile

/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/xarray/core/duck_array_ops.py:675: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
  return push(array, n, axis)

<xarray.Dataset>
Dimensions:      (cycle: 6, σ: 60, bounds: 2)
Coordinates:
  * cycle        (cycle) int64 0 1 2 3 4 5
  * σ            (σ) float32 34.15 34.16 34.17 34.18 ... 37.2 37.2 37.2 37.2
    z_σ          (cycle, σ) float64 66.9 67.23 67.65 ... 5.036e+03 5.166e+03
    z_σ_bounds   (cycle, σ, bounds) float64 66.73 67.06 ... 5.101e+03 5.232e+03
Dimensions without coordinates: bounds
Data variables:
    KAPPA_ISOP   (cycle, σ) float64 [m²/s] 2.488e+03 2.486e+03 ... 508.3 464.9
    RediVar      (cycle, σ) float64 [K²/s] 2.956e-09 2.951e-09 ... 8.636e-13
    TEMP         (cycle, σ) float64 [K] 295.8 295.8 295.7 ... 274.8 274.8 274.8
    delT2        (cycle, σ) float64 [K²/m²] 1.157e-12 1.156e-12 ... 1.605e-15
    delT2_plane  (cycle, σ) float64 [K²/m²] 1.157e-12 1.156e-12 ... 1.605e-15
Attributes:
    Conventions:       CF-1.0; http://www.cgd.ucar.edu/cms/eaton/netcdf/CF-cu...
    calendar:          All years have exactly  365 days.
    cell_methods:      cell_methods = time: mean ==> the variable values are ...
    contents:          Diagnostic and Prognostic Variables
    history:           none
    model_doi_url:     https://doi.org/10.5065/D67H1H0V
    revision:          $Id$
    source:            CCSM POP2, the CCSM Ocean Component
    start_time:        This dataset was created on 2019-06-20 at 13:53:17.7
    time_period_freq:  month_1
    title:             g.e21.GOMIPECOIAF_JRA.TL319_g17.CMIP6-omip2.001

POP 1° spinup#

%autoreload

ds_ = xr.load_dataset(
    "../datasets/pop-1deg-redi-var-natre-cycle-0.zarr/",
    engine="zarr",
)
ds_.z_σ.attrs["units"] = "cm"
ds_.σ.attrs.update({"axis": "Z", "positive": "down"})
ds_.TEMP.attrs["standard_name"] = "sea_water_potential_temperature"

month1 = ds_.isel(yearmonth=0).pint.quantify().map(dcpy.util.to_base_units).squeeze()
month1_profile = ed.pop.calc_mean_redivar_profile(month1)
month1_profile = month1_profile.where(month1_profile.z_σ.notnull(), drop=True).isel(
    σ=slice(1, None)
)

pop_1deg_spinup = (
    ds_.squeeze("cycle")
    .coarsen(yearmonth=120, boundary="trim")
    .mean()
    .rename({"yearmonth": "decade"})
    .pint.quantify()
    .map(dcpy.util.to_base_units)
)
pop_1deg_spinup_profile = ed.pop.calc_mean_redivar_profile(pop_1deg_spinup)
pop_1deg_spinup_profile

/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/xarray/core/duck_array_ops.py:675: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
  return push(array, n, axis)
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/xarray/core/duck_array_ops.py:675: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
  return push(array, n, axis)

<xarray.Dataset>
Dimensions:      (decade: 6, σ: 60, bounds: 2)
Coordinates:
    cycle        int64 0
    month        (decade) float64 6.5 6.5 6.5 6.5 6.5 6.5
    time         (decade) object 0005-12-16 00:00:00 ... 0055-12-16 00:00:00
    year         (decade) float64 5.5 15.5 25.5 35.5 45.5 55.5
  * σ            (σ) float32 34.15 34.15 34.17 34.18 ... 37.2 37.2 37.2 37.2
    z_σ          (decade, σ) float64 66.59 66.93 67.3 ... 4.931e+03 4.931e+03
    z_σ_bounds   (decade, σ, bounds) float64 66.42 66.76 ... 4.931e+03 4.931e+03
Dimensions without coordinates: decade, bounds
Data variables:
    KAPPA_ISOP   (decade, σ) float64 [m²/s] 2.512e+03 2.518e+03 ... nan nan
    RediVar      (decade, σ) float64 [K²/s] 3.03e-09 3.025e-09 ... nan nan
    TEMP         (decade, σ) float64 [K] 295.8 295.8 295.7 295.7 ... nan nan nan
    delT2        (decade, σ) float64 [K²/m²] 1.186e-12 1.183e-12 ... nan nan
    delT2_plane  (decade, σ) float64 4.255e-13 4.351e-13 4.502e-13 ... 0.0 0.0
Attributes:
    Conventions:       CF-1.0; http://www.cgd.ucar.edu/cms/eaton/netcdf/CF-cu...
    calendar:          All years have exactly  365 days.
    cell_methods:      cell_methods = time: mean ==> the variable values are ...
    contents:          Diagnostic and Prognostic Variables
    history:           none
    model_doi_url:     https://doi.org/10.5065/D67H1H0V
    revision:          $Id$
    source:            CCSM POP2, the CCSM Ocean Component
    start_time:        This dataset was created on 2019-06-20 at 13:53:17.7
    time_period_freq:  month_1
    title:             g.e21.GOMIPECOIAF_JRA.TL319_g17.CMIP6-omip2.001

Make DataTree#

pop_natre = DataTree()
pop_natre["1deg/cycle"] = DataTree(pop_1deg_natre_profile)
pop_natre["hires/clim"] = DataTree(
    xr.merge([pop_hires_natre, pop_hires_clim_profile], join="exact")
)
pop_natre["1deg/decadal"] = DataTree(pop_1deg_spinup_profile)
pop_natre["1deg/month1"] = DataTree(month1_profile)
pop_natre

<xarray.DatasetView>
Dimensions:  ()
Data variables:
    *empty*

Cole/Groeskamp#

This uses a crappier gradient.

Hide code cell source 
cole = ed.read_cole()

# read, take mean in region; swap to pressure coordinate
argograd = (
    xr.open_dataset("../datasets/cole-clim-gradient-natre.nc")
    .reset_coords("pres")
    .sel(lat=slice(24, 28), lon=slice(360 - 30, 360 - 25))
    .cf.guess_coord_axis()
)
argograd.lon.attrs["axis"] = "X"
argograd["lon"] = argograd.lon.cf.guess_coord_axis()
counts = argograd.CT.count(["lat", "lon"])
argo_delT2_plane = ed.plane_fit_gradient(argograd.CT).assign_coords(
    {"pres": argograd.pres.mean(["lon", "lat"])}
)

cole_natre = cole[["diffusivity", "density_mean_depth", "depth_mean_sig"]].interp(
    lat=argograd.lat.data, lon=argograd.lon.data
)

# move from sigma to depth space
argo_dTiso = xgcm.transform.linear_interpolation(
    argograd.dTiso,
    argograd.pres,
    cole.pres,
    "sigma",
    "sigma",
    "pres",
)

cole_natre["RediVar"] = cole_natre.diffusivity * argo_dTiso**2
cole_natre["RediVar"].attrs = {
    "long_name": "$K_e^{cole} |∇T_{iso}^{argo}|²$",
    "units": "°C²/s",
}
cole_natre = cole_natre.mean(["lat", "lon"])

groeskamp = (
    xr.open_dataset("~/datasets/groeskamp2020/groeskamp2020.nc")
    .sel(lat=slice(24, 28), lon=slice(360 - 30, 360 - 25))
    .rename({"depth": "pres"})
)
groeskamp["dTiso"] = argo_dTiso.interp(pres=groeskamp.pres)
groeskamp["eddy_var_0"] = groeskamp["Ke_0"] * groeskamp.dTiso**2
groeskamp["eddy_var"] = groeskamp["Ke"] * groeskamp.dTiso**2
groeskamp = groeskamp.mean(["lat", "lon"])

/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/xarray/core/indexing.py:661: RuntimeWarning: deallocating CachingFileManager(<class 'netCDF4._netCDF4.Dataset'>, '/Users/dcherian/datasets/groeskamp2020/groeskamp2020.nc', mode='r', kwargs={'clobber': True, 'diskless': False, 'persist': False, 'format': 'NETCDF4'}, manager_id='641c78cd-86a0-45ef-90b8-0a944d6638e8'), but file is not already closed. This may indicate a bug.
  self.array = NumpyIndexingAdapter(np.asarray(self.array))

a05 = xr.load_dataset("../datasets/a05-2015.nc")
a05

<xarray.Dataset>
Dimensions:           (bound: 2, gamma_n: 13, cast: 46)
Coordinates: (12/13)
  * bound             (bound) object 'lower' 'upper'
  * gamma_n           (gamma_n) float64 26.6 26.88 27.11 ... 27.93 27.96 27.98
    num_obs           (gamma_n) int64 3120 3281 3326 3256 ... 3549 3740 3463
    pres              (gamma_n) float64 233.1 383.6 ... 1.976e+03 2.149e+03
    gamma_n_bounds    (bound, gamma_n) float64 26.45 26.75 27.0 ... 27.97 27.99
    station           (cast) int64 108 109 110 111 112 ... 126 127 128 129 130
    ...                ...
    cleaner           <U2 'GE'
    time              (cast) datetime64[ns] 2016-01-09T06:22:19.077246976 ......
    lon               (cast) float64 325.1 325.7 326.4 ... 337.7 338.4 339.0
    lat               (cast) float64 24.5 24.5 24.5 24.5 ... 24.92 25.13 25.33
    CTD_chipod        float64 9.969e+36
  * cast              (cast) int64 0 1 2 3 4 5 6 7 8 ... 38 39 40 41 42 43 44 45
Data variables: (12/59)
    chi               (gamma_n) float64 1.788e-08 7.233e-09 ... 2.8e-10
    eps               (gamma_n) float64 2.916e-09 1.259e-09 ... 2.938e-10
    KtTz              (gamma_n) float64 2.282e-07 2.495e-07 ... 6.227e-08
    KtTz~             (gamma_n) float64 3.673e-07 1.921e-07 ... 6.11e-08
    chib2             (gamma_n) float64 8.941e-09 3.617e-09 ... 1.4e-10
    hm                (gamma_n) float64 148.4 149.3 151.3 ... 165.5 173.1 162.3
    ...                ...
    δKtTzTz           (gamma_n) float64 2.53e-09 2.664e-09 ... 5.358e-11
    δKtTz~Tz          (gamma_n) float64 1.783e-09 9.228e-10 ... 7.829e-11
    δwTTz             (gamma_n) float64 nan nan nan nan nan ... nan nan nan nan
    δresidual         (gamma_n) float64 -1.123e-09 -3.266e-10 ... -1.414e-11
    δresidual_chi     (gamma_n) float64 4.542e-10 1.32e-10 ... 2.716e-11
    δpres             (gamma_n) float64 74.2 74.67 75.67 ... 82.74 86.57 81.13
Attributes: (12/15)
    cruise_name:       A05
    featureType:       trajectoryProfile
    Conventions:       CF-1.8
    title:             Turbulence quantities measured by CTD-chipods mounted ...
    summary:           Turbulence data across the ocean basin from 10 m to th...
    processing_level:  3: Variables mapped on uniform space-time grid from hi...
    ...                ...
    creator_name:      Jonathan Nash / Aurelie Moulin 
    institution:       Ocean Mixing Group, College of Earth, Ocean, and Atmos...
    project:           A05, GO-SHIP (The Global Ocean Ship-Based Hydrographic...
    reference1:        J.D. Nash et al., Ocean mixing measured by fast-respon...
    comment:           Expocode: 740H20200119
    commit:            eddydiff: b55f9526739dd63a768df2f816bc8327bcc0d247 |  ...

MIMOC#

mimoc = (
    dcpy.oceans.read_mimoc(globstr="MIMOC_SIGMA0_*")
    .coarsen(latitude=2, longitude=2, boundary="trim")
    .mean()
)
mimoc.coords["lat_face"] = (
    "lat_face",
    mimoc.latitude.data - 0.25,
    {"axis": "Y", "c_grid_axis_shift": -0.5},
)
mimoc.coords["lon_face"] = (
    "lon_face",
    mimoc.longitude.data - 0.25,
    {"axis": "X", "c_grid_axis_shift": -0.5},
)

dx = gsw.distance(mimoc.lon_face.data, np.zeros_like(mimoc.lon_face.data))
dx = np.append(dx, dx[-1])
mimoc.coords["dx"] = ("longitude", dx, {"units": "m"})
mimoc.coords["dx"] = mimoc.dx * np.cos(np.pi / 180 * mimoc.latitude)
dy = gsw.distance(np.zeros_like(mimoc.latitude.data), mimoc.lat_face.data)
dy = np.append(dy, dy[-1])
mimoc.coords["dy"] = ("latitude", dy, {"units": "m"})
mimoc.coords["X"] = mimoc.dx.cf.cumsum("X")
mimoc.coords["Y"] = mimoc.dy.cf.cumsum("Y")

grid = xgcm.Grid(mimoc, metrics={"X": ["dx"], "Y": ["dy"]}, periodic=["X", "Y"])
dTdx = grid.derivative(mimoc.CT, "X")
dTdy = grid.derivative(mimoc.CT, "Y")
mimoc["delT2"] = grid.interp(dTdx, "X") ** 2 + grid.interp(dTdy, "Y") ** 2
mimoc.delT2.attrs = {"long_name": "$|∇T|^2$", "units": "C/m"}

mimoc_natre = mimoc.sel(
    latitude=slice(25, 30), longitude=slice(360 - 30, 360 - 25)
).load()
mimoc_natre_profile = (
    mimoc_natre.reset_coords("pressure")
    .cf.mean(["X", "Y", "time"])
    .set_coords("pressure")
)
mimoc_natre.sigma0.attrs["axis"] = "Z"
mimoc

/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/xgcm/grid.py:1645: UserWarning: Metric at ('time', 'lon_face', 'latitude', 'sigma0') being interpolated from metrics at dimensions ('longitude', 'latitude'). Boundary value set to 'extend'.
  warnings.warn(
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/xgcm/grid.py:1645: UserWarning: Metric at ('time', 'longitude', 'lat_face', 'sigma0') being interpolated from metrics at dimensions ('latitude',). Boundary value set to 'extend'.
  warnings.warn(

<xarray.Dataset>
Dimensions:        (time: 12, longitude: 360, latitude: 170, sigma0: 75,
                    lat_face: 170, lon_face: 360)
Coordinates:
  * latitude       (latitude) float32 -79.75 -78.75 -77.75 ... 87.25 88.25 89.25
  * longitude      (longitude) float32 0.25 1.25 2.25 3.25 ... 357.2 358.2 359.2
  * sigma0         (sigma0) float32 22.4 22.64 22.88 23.11 ... 28.14 28.15 28.24
  * time           (time) datetime64[ns] 2014-01-15 2014-02-15 ... 2014-12-15
  * lat_face       (lat_face) float32 -80.0 -79.0 -78.0 -77.0 ... 87.0 88.0 89.0
  * lon_face       (lon_face) float32 0.0 1.0 2.0 3.0 ... 357.0 358.0 359.0
    dx             (longitude, latitude) float64 1.979e+04 ... 1.455e+03
    dy             (latitude) float64 1.112e+05 1.112e+05 ... 1.112e+05
    X              (longitude, latitude) float64 1.979e+04 ... 5.24e+05
    Y              (latitude) float64 1.112e+05 2.224e+05 ... 1.879e+07 1.89e+07
Data variables:
    SUMMED_WEIGHT  (time, longitude, latitude, sigma0) float32 dask.array<chunksize=(1, 360, 170, 75), meta=np.ndarray>
    YEAR_OF_DATA   (time, longitude, latitude, sigma0) float32 dask.array<chunksize=(1, 360, 170, 75), meta=np.ndarray>
    pressure       (time, longitude, latitude, sigma0) float32 dask.array<chunksize=(1, 360, 170, 75), meta=np.ndarray>
    SA             (time, longitude, latitude, sigma0) float32 dask.array<chunksize=(1, 360, 170, 75), meta=np.ndarray>
    CT             (time, longitude, latitude, sigma0) float32 dask.array<chunksize=(1, 360, 170, 75), meta=np.ndarray>
    delT2          (time, longitude, latitude, sigma0) float64 dask.array<chunksize=(1, 360, 170, 75), meta=np.ndarray>
mimoc_natre.CT.isel(sigma0=24, time=0).cf.plot()

<matplotlib.collections.QuadMesh at 0x178bd8e20>

mimoc_natre_profile["delT2_plane"] = ed.plane_fit_gradient(
    mimoc_natre.CT.mean("time"), debug=True
)
argo_delT2_plane = ed.plane_fit_gradient(argoclim.CT, debug=True).assign_coords(
    {"pres": argoclim.pres.mean(["lon", "lat"])}
)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[23], line 4
      1 mimoc_natre_profile["delT2_plane"] = ed.plane_fit_gradient(
      2     mimoc_natre.CT.mean("time"), debug=True
      3 )
----> 4 argo_delT2_plane = ed.plane_fit_gradient(argoclim.CT, debug=True).assign_coords(
      5     {"pres": argoclim.pres.mean(["lon", "lat"])}
      6 )

NameError: name 'argoclim' is not defined

AVISO#

aviso = xr.open_zarr("~/datasets/aviso_monthly.zarr", consolidated=True)
aviso["longitude"] = aviso.longitude - 360
aviso = aviso.sel(latitude=slice(10, 50), longitude=slice(-40, 5))
aviso

Variance production vs dissipation#

  1. Agreement between POP “mesoscale variance dissipation” and microstructure estimate of “variance received from mesoscale” is great! (Purple dashed steps agree well with purple solid line)

  2. POP’s grid-scale temp variance dissipation is really dominated by the horizontal biharmonic diffusion, but fundamentally the amount of mesoscale variance destroyed by grid-scale diffusion matches that in the microstructure measurements.

  3. One place of disagreement is between 400m and 700m where the POP estimate is really small compared to my difference. FP2005 also see this using the mooring based Ke estimate. At these levels their two terms also match a lot closer than mine, so there’s an error somewhere. It seems like my \(⟨χ⟩/2\) is a little bit too high which happened after I switched to the bootstrapping based estimates. I’ve made some improvements

    1. Sorting the dataset, changes the mean density field used to pick the bins.

    2. Calculating \((∂T/∂z)\_m\) is a little wonky. I now use a different procedure. In each O(100m) \(γ_n\) bin, I groupby finer-resolution \(γ_n\) bins (21 of them), calculate the mean profile, then fit a straight line.

  4. Cole/Groeskamp \(K_e\) estimate based variance dissipation compares much better now.

    1. ~Comparing to variance dissipated by parameterizations for \(K_e\) (\(K_e |∇_hT|²\)) doesn’t work well. I wonder if I’m doing something wrong with my \(∇_hT\) estimate~

    2. I updated my |∇T| estimation code; now estimate along-isopycnal gradients in σ space instead of just a lateral gradient

    3. While the current algo still doesn’t match Silvia’s ∇S (but much better than earlier), it seems to improve agreement (panel c)

    4. Still a lot of disagreement below 1200m when using published \(K_e\) estimates.

Hide code cell source 
    plt.rcParams["figure.dpi"] = 140
plt.rcParams["figure.facecolor"] = "none"  # (0.91, 0.91, 0.91)
plt.rcParams["axes.facecolor"] = (0.99,) * 3
plt.rcParams.update(
    {
        "axes.spines.top": True,
        "axes.spines.bottom": True,
        "axes.spines.left": True,
        "axes.spines.right": True,
        "axes.edgecolor": (0.31,) * 3,
        "axes.grid.which": "both",
        "grid.linewidth": 0.5,
        "grid.alpha": 0.3,
        "grid.color": (0,) * 3,
    }
)


micro.pres.attrs["bounds"] = "pres_err"
micro["residual"] = np.abs(micro.residual)
micro["residual_err"] = np.abs(micro.residual_err)
hires = pop_natre["hires/clim"].drop_vars(["z_σ", "z_σ_bounds"])

mesomicrocolor = "darkorchid"
error_kwargs = dict(error="x", ls="none", marker=".")


def plot_pop_natre_meso_micro(ax):
    (-1 * hires["DISS"]).cf.plot(
        ax=ax, lw=1.5, color=mesomicrocolor, label="POP meso $→$ micro"
    )
    dcpy.plots.errorbar(
        micro,
        x="residual",
        y="pres",
        color=mesomicrocolor,
        label="NATRE: $⟨χ⟩/2 - K_ρ^m ∂_zθ_m^2$",
        # lw=2.5,
        ax=ax,
        **error_kwargs,
    )

    # dcpy.plots.fill_between_bounds(
    #    micro,
    #    "residual",
    #    y="pres",
    #    color=mesomicrocolor,
    #    label="NATRE meso $→$ micro",
    #    ax=ax,
    # )


# f, axx = plt.subplots(1, 4, )
f, ax = plt.subplot_mosaic(
    [["natre", "pop-hires", "params"]],
    sharey=False,
    sharex=True,
    constrained_layout=False,
)

### Mesoscale budget
hires["BC"].cf.plot(ax=ax["pop-hires"], label="POP h. stir $⟨u_e^h θ_e⟩.∇_hθ_m$")
np.abs(-1 * hires["PKC"]).cf.plot(
    ax=ax["pop-hires"], label="POP v. stir $-⟨w_eθ_e⟩.∂_zθ_m$"
)
# (-1 * pop_natre.VMIX).cf.plot(color="sienna", ax=ax[0], label="POP VMIX")
# (-1 * pop_natre.HDIFF).cf.plot(ax=ax[0], label="POP HDIFF", marker="x", ls="none")


#################################
######## Microscale budget: NATRE
#################################
dcpy.plots.fill_between_bounds(
    micro,
    "KρTz2",
    y="pres",
    color="k",
    label="FP2005: $K_ρ^m ∂_zθ_m^2$",
    ax=ax["natre"],
)
# dcpy.plots.fill_between_bounds(
#    micro,
#    "KtTz~Tz",
#    y="pres",
#    color="C1",
#    label="NATRE $⟨χ/2/∂_z\widetilde{T}⟩ ∂_zθ_m$",
#    ax=ax[1],
# )
dcpy.plots.fill_between_bounds(
    micro,
    "chib2",
    y="pres",
    color="C0",
    label="FP2005: $⟨χ⟩/2$",
    ax=ax["natre"],
)


# (argo_fine_depth.χ.sel(criteria="whalen") / 2).cf.plot.step(
#    # marker="o",
#    color="green",
#    ax=ax[1],
#    y="pressure_bins",
#    label="$⟨χ⟩^{argo}_{Rω=3}/2$",
# )

#################################
########## Microscale budget: A05
#################################
if "a05" in ax:
    dcpy.plots.fill_between_bounds(
        a05_binned,
        "KtTz~Tz",
        y="pres",
        color="k",
        label=r"$⟨χ/2/∂_z\widetilde{T}⟩ ∂_zθ_m$",
        ax=ax["a05"],
    )
    dcpy.plots.fill_between_bounds(
        a05_binned,
        "chib2",
        y="pres",
        color="C0",
        label="$⟨χ⟩/2$",
        ax=ax["a05"],
    )
    dcpy.plots.errorbar(
        a05_binned,
        x="residual_chi",
        y="pres",
        color=mesomicrocolor,
        label="A05 meso → micro",
        ax=ax["a05"],
        **error_kwargs,
    )

#################################
### Parameterizations
#################################
cole_natre.RediVar.plot(
    y="pres", label="$K_e^{cole} |∇_hT^{cole}|²$", color="k", ls="--", ax=ax["params"]
)

# dcpy.plots.errorbar(
#    micro,
#    x="residual",
#    y="pres",
#    # color=mesomicrocolor,
#    label="NATRE $⟨χ⟩/2 - ⟨χ/2/∂_z\widetilde{T}⟩ ∂_zθ_m$",
#    ax=ax[3],
#    **error_kwargs,
# )
if "a05" in ax:
    dcpy.plots.errorbar(
        a05_binned,
        x="residual_chi",
        y="pres",
        # color=mesomicrocolor,
        label=r"A05 $⟨χ⟩/2 - ⟨χ/2/∂_z\widetilde{T}⟩ ∂_zθ_m$",
        ax=ax["params"],
        **error_kwargs,
    )

# groeskamp.eddy_var_0.plot(
#    y="pres",
#    label="$K_e^{G2020unsupp} |∇_hT^{cole}|²$",
#    color="limegreen",
#    ls="-.",
#    ax=ax[2],
# )
groeskamp.eddy_var.plot(
    y="pres",
    label="$K_e^{G2020} |∇_hT^{cole}|²$",
    color="limegreen",
    ls="--",
    ax=ax["params"],
)

#################################
### Cleanup
#################################
for axx in ax.values():
    plot_pop_natre_meso_micro(axx)
    axx.set_xscale("log")
    axx.set_xlabel("")
    axx.set_ylim([1850, 350])
    axx.set_xlim([1e-11, 1e-8])
    # axx.set_xticks([1e-11, 1e-10, 1e-9, 1e-8])
    axx.xaxis.set_minor_locator(
        mpl.ticker.LogLocator(base=10, subs="all", numticks=120)
    )
    #  axx.set_ylim([1900, 180])
    axx.legend(loc="upper center", bbox_to_anchor=(0.5, -0.15))
    axx.grid(True, which="both", lw=0.5)

    patch = mpl.patches.Rectangle(
        (0, 800), 1, 700, transform=axx.get_yaxis_transform(), alpha=0.1, zorder=1
    )
    axx.add_patch(patch)
# ax[1].set_xlabel("Variance production or dissipation [°C²/s]")

# f.suptitle("NATRE region")
ax["pop-hires"].set_title("Mesoscale budget: POP2 1/10°", fontsize="medium")
ax["natre"].set_title("Microscale budget: NATRE", fontsize="medium")
if "a05" in ax:
    ax["a05"].set_title("Microscale budget: A05", fontsize="medium")
ax["params"].set_title("$K_e$ Parameterizations", fontsize="medium")


dcpy.plots.label_subplots(ax.values(), x=0.03, y=0.95)
dcpy.plots.clean_axes(list(ax.values()))
f.set_size_inches((8, 5))
f.savefig("../images/natre-meso-micro-param.pdf", bbox_inches="tight")
f.savefig("../images/natre-meso-micro-param.png", bbox_inches="tight")
_images/96c7b6109e718715519e4c7dc17c55798ed7a936eb7ea860db2c2affb73239e4.png
Hide code cell source 
micro.pres.attrs["bounds"] = "pres_err"

mesomicrocolor = "darkorchid"
error_kwargs = dict(error="x", ls="none", marker=".")


def plot_pop_natre_meso_micro(ax):
    (-1 * pop_natre.DISS).cf.plot(
        ax=ax, lw=1.5, color=mesomicrocolor, label="POP meso $→$ micro"
    )

    # dcpy.plots.fill_between_bounds(
    #    micro,
    #    "residual",
    #    y="pres",
    #    color=mesomicrocolor,
    #    label="NATRE meso $→$ micro",
    #    ax=ax,
    # )


f, ax = plt.subplots(1, 3, sharey=False, sharex=True, constrained_layout=False)

### Mesoscale budget
# pop_natre.BC.cf.plot(ax=ax[0], label="POP h. stir $⟨u_e^h θ_e⟩.∇_hθ_m$")
# np.abs(-1 * pop_natre.PKC).cf.plot(ax=ax[0], label="POP v. stir $-⟨w_eθ_e⟩.∂_zθ_m$")
# (-1 * pop_natre.VMIX).cf.plot(color="sienna", ax=ax[0], label="POP VMIX")
# (-1 * pop_natre.HDIFF).cf.plot(ax=ax[0], label="POP HDIFF", marker="x", ls="none")


#################################
######## Microscale budget: NATRE
#################################
plt.sca(ax[0])

dcpy.plots.fill_between_bounds(
    micro,
    "KρTz2",
    y="pres",
    color="k",
    label="FP2005: $χ_m = K_ρ^m ∂_zθ_m^2$",
    ax=ax[0],
)
# dcpy.plots.fill_between_bounds(
#    micro,
#    "KtTz~Tz",
#    y="pres",
#    color="C1",
#    label="NATRE $⟨χ/2/∂_z\widetilde{T}⟩ ∂_zθ_m$",
#    ax=ax[0],
# )
dcpy.plots.fill_between_bounds(
    micro,
    "chib2",
    y="pres",
    color="C0",
    label="FP2005: $⟨χ⟩/2$",
    ax=ax[0],
)

dcpy.plots.errorbar(
    micro,
    x="residual",
    y="pres",
    color=mesomicrocolor,
    label="NATRE $χ_e$",
    ax=ax[0],
    **error_kwargs,
)

# (argo_fine_depth.χ.sel(criteria="whalen") / 2).cf.plot.step(
#    # marker="o",
#    color="green",
#    ax=ax[1],
#    y="pressure_bins",
#    label="$⟨χ⟩^{argo}_{Rω=3}/2$",
# )

#################################
########## Microscale budget: A05
#################################
dcpy.plots.fill_between_bounds(
    a05_binned,
    "KtTz~Tz",
    y="pres",
    color="k",
    label="$χ_m$",
    ax=ax[1],
)
dcpy.plots.fill_between_bounds(
    a05_binned,
    "chib2",
    y="pres",
    color="C0",
    label="$⟨χ⟩/2$",
    ax=ax[1],
)
dcpy.plots.errorbar(
    a05_binned,
    x="residual_chi",
    y="pres",
    color=mesomicrocolor,
    label="A05 $χ_e$",
    ax=ax[1],
    **error_kwargs,
)

#################################
### Parameterizations
#################################
cole_natre.RediVar.plot(
    y="pres", label="$K_e^{cole} |∇_hT^{cole}|²$", color="k", ls="--", ax=ax[2]
)
dcpy.plots.errorbar(
    micro,
    x="residual",
    y="pres",
    # color=mesomicrocolor,
    label="NATRE: $χ_e$",
    lw=2.5,
    ax=ax[2],
    **error_kwargs,
)
# dcpy.plots.errorbar(
#    micro,
#    x="residual",
#    y="pres",
#    # color=mesomicrocolor,
#    label="NATRE $⟨χ⟩/2 - ⟨χ/2/∂_z\widetilde{T}⟩ ∂_zθ_m$",
#    ax=ax[3],
#    **error_kwargs,
# )
dcpy.plots.errorbar(
    a05_binned,
    x="residual_chi",
    y="pres",
    # color=mesomicrocolor,
    label="A05 $χ_e$",
    ax=ax[2],
    **error_kwargs,
)

# groeskamp.eddy_var_0.plot(
#    y="pres",
#    label="$K_e^{G2020unsupp} |∇_hT^{cole}|²$",
#    color="limegreen",
#    ls="-.",
#    ax=ax[2],
# )
groeskamp.eddy_var.plot(
    y="pres", label="$K_e^{G2020} |∇_hT^{cole}|²$", color="limegreen", ls="--", ax=ax[2]
)

#################################
### Cleanup
#################################
for axx in ax:
    plot_pop_natre_meso_micro(axx)
    axx.set_xscale("log")
    axx.set_xlabel("")
    axx.set_ylim([1850, 350])
    axx.set_xlim([1e-11, 1e-8])
    # axx.set_xticks([1e-11, 1e-10, 1e-9, 1e-8])
    axx.xaxis.set_minor_locator(
        mpl.ticker.LogLocator(base=10, subs="all", numticks=120)
    )
    #  axx.set_ylim([1900, 180])
    axx.legend(loc="upper center", bbox_to_anchor=(0.5, -0.15))
    axx.grid(True, which="both", lw=0.5)
# ax[1].set_xlabel("Variance production or dissipation [°C²/s]")

f.suptitle("NATRE region")
# ax[0].set_title("(a) Mesoscale budget", fontsize="medium")
ax[0].set_title("(a) Microscale budget: NATRE", fontsize="medium")
ax[1].set_title("(b) Microscale budget: A05", fontsize="medium")
ax[2].set_title("(c) Parameterizations", fontsize="medium")

dcpy.plots.clean_axes(ax)
f.set_size_inches((8, 5))
f.savefig("../images/natre-3-panel.png", bbox_inches="tight")

Interpretation#

Our analytical framework follows that of Ferrari & Polzin (2005) and Garrett (1992). Begin with a triple decomposition of the temperature fields

(2)#\[\begin{equation} T = T_m + T_e + T_t \end{equation}\]

where the subscripts \(m\), \(e\) and \(t\) are loosely identified with the “mean”, “mesoscale eddy”, and “turbulence” scales.

Define the “large-scale averaging operator” \(⟨⟩\) that filters out the “eddy” and “turbulence” so that \(⟨T⟩ = T_m\) and an intermediate scale averaging operator \(\widetilde{\phantom{...}}\) that filters out the “turbulence” so that \(\widetilde{T} = T_m + T_e\).

Ferrari & Polzin (2005) follow Garrett (1992) to derive the variance budget equations for \(⟨T_e^2⟩\) and \(⟨T_t^2⟩\). Assuming stationarity, and homogeneity, and /a priori/ ignoring the triple correlation term, we can write the approximate equations

(3)#\[\begin{align} ⟨u_eT_e⟩ ⋅ ∇T_m - ⟨\widetilde{u_t T_t} ⋅ ∇T_e⟩ &\approx 0 \\ ⟨\widetilde{u_tT_t} ⋅ ∇ (T_m + T_e)⟩ &\approx - \frac 12 ⟨\widetilde{χ}⟩; \qquad χ = κ_T |∇T_t|^2 \end{align}\]

Here \(χ\) representes ths instantaneous rate of dissipation of termpeature variance and \(κ_T\) is the molecular diffusiovity of temperature.

The first equation, for mesoscale variance \(⟨T_e^2⟩\), states that the mesoscale scale field generates variance by the stirring of the mean and this variance cascades down to the microscale turbulence through the “scale transoframtion term” \(⟨\widetilde{u_t T_t} ⋅ ∇T_e⟩\). The second equation, for microscale variance \(⟨T_t^2⟩\), states that the microscale turbulence stirs the combination of the mean and eddy fields \((T_m + T_e)\) to generate variance that is eventually dissipated at the moelcular scale at the rate \(χ\). The scale transformation term \(⟨\widetilde{u_t T_t} ⋅ ∇T_e⟩ \)has opposite sign in the two equations indicating that it is a link between the mesoscale and microscale. Adding the two equations yields the approximate blaance,

(4)#\[\begin{equation} ⟨\widetilde{u_tT_t}⟩ ⋅ ∇ T_m + ⟨\widetilde{u_e T_e}⟩ ⋅ ∇ T_m &\approx - \frac 12 ⟨χ⟩ \end{equation}\]

that is the variance dissipated at the molecular scale is approximately generated by the stirring of the mean field \(T_m\) by the mesoscale and microscale turbulence.

A large number of approximations are required to get to this point, but this framework has value in qualitatively describing the nature of variance cascades in the ocean \citep{Garrett1992,Davis1994,Ferrari2005,NaveriaGarabato2016,Springys2020}.

(5)#\[\begin{align} ⟨\widetilde{u_tT_t}⟩ ⋅ ∇ T_m + ⟨\widetilde{u_t T_t} ⋅ ∇ T_e⟩ &\approx - \frac 12 ⟨χ⟩ \\ ⟨\widetilde{u_tT_t}⟩ ⋅ ∇ T_m + ⟨u^h_eT_e⟩ ⋅ ∇^hT_m &\approx - \frac 12 ⟨χ⟩ \end{align}\]
(6)#\[\begin{align} ⟨\widetilde{u_tT_t}⟩ &\approx ⟨w_t T_t⟩ \approx ⟨w_tb_t⟩ = \frac {Γ ⟨ε⟩}{N_m^2} \\ ⟨\widetilde{u_tT_t}⟩ &\approx ⟨\widetilde{w_t T_t}⟩ \approx -\frac 12 \left\langle\frac{χ}{∂_z\widetilde T}\right\rangle \end{align}\]

Spicy region#

We can interpret the red mesoscale dissipation term in (b) as the amount of variance transferred to the microscale for eventual dissipation. This peaks in regions of high mean spice between 1000m and 1500m. The balance in spicy regions is that horizontal mesoscale stirring converts MPE → EPE (compensated filaments or T-S variance) that is dissipated by the microscale. The approximate budget is (note vertical stirring gets ignored, and \(u^h\) represents horizontal velocity vector)

(7)#\[\begin{align} ⟨u^h_eT_e⟩ ⋅ ∇_hT_m - ⟨\widetilde{u_t T_t} ⋅ ∇T_e⟩ &\approx 0 \\ ⟨\widetilde{u_t T_t} ⋅ ∇ T_e⟩ &\approx - \frac 12 ⟨χ⟩ \end{align}\]

  1. Yiming sees the first equation as “lateral stirring” balances “mesoscale dissipation”, and

  2. FP2005 combine the two equations to conclude that “mesoscale stirring of mean variance” balances “total microscale dissipation”

Weak spice#

In the upper water column that is not spicy, we instead see horizontal stirring MPE→ EPE which is then routed to EKE by “vertical stirring”. Here the scale transformation term is 0 so the budget is

(8)#\[\begin{align} ⟨u_eT_e⟩ ⋅ ∇T_m &\approx 0 \\ ⟨u_tT_t⟩ ⋅ ∇ T_m &\approx - \frac 12 ⟨χ⟩ \end{align}\]

  1. Yiming sees the first equation as “low dissipation” or “lateral stirring balances vertical stirring”;

  2. FP2005 sees the second equation as “microscale stirring of the mean” balances “total microscale dissipation”

All microscale budget terms#

Microscale budget:

  • Ferrari & Polzin use \(⟨u_t θ_t⟩.∇θ_m = -Γ ⟨ε⟩/N_m^2 ∂_zθ_m^2\), we are instead doing \(-⟨K_T θ_z⟩ ∂_zθ_m\).

  • Our estimate mostly agrees with theirs (crosses versus solid-line-steps).

    1. Filling in NaNs in the temperature profiles was important. This means more estimates are used in the averages

    2. Masked where both \(χ\) and \(ε\) are both available.

    3. Masked so that \(|θ_z|\) > 1e-3 (really important, otherwise our estimate is too high)

Variance produced by parameterizations for eddy diffusivity:

  • The comparison with variance production using Cole et al diffusivity and gradients calculated in a similar way disagree significantly between 400 and 1000m. Below that there’s not much data but it seems to agree.

  • The Groeskamp et al (2020) estimate [G2020] is slightly better but not by much. All predict a maximum in variance production near 800m

⟨χ⟩ estimate using finescale parameterizations:

  • The finescale estimates are interesting, argo with Rω=3 and natre with Rω=7 work well but why?

  • There is some interesting discrepancy here. This choice of \(Rω\) makes for a factor of 2 difference which would hide the mesoscale→microscale transformation.

f, ax = plt.subplots(1, 1)  # , constrained_layout=True)

(micro.chi / 2).cf.plot.step(
    y="pres",
    color="r",
    label="FP: $⟨χ⟩/2$",
    lw=2,
    xlim=(1e-11, 1e-8),
    ylim=(1800, 200),
)
(micro.Krho_m * micro.dTdz_m**2).cf.plot.step(
    color="k", lw=2, y="pres", label="FP: $Γ ⟨ε⟩/N_m^2  ∂_zθ_m^2$"
)

cole_natre.RediVar.plot(y="pres", label="$K_e^{cole} |∇T^{cole}|²$", color="k", ls="--")

groeskamp.eddy_var_0.plot(
    y="depth", label="$K_e^{G2020unsupp} |∇T^{cole}|²$", color="teal", ls="-."
)
groeskamp.eddy_var.plot(
    y="pres", label="$K_e^{G2020} |∇T^{cole}|²$", color="teal", ls="--"
)

ed.sections.plot_var_prod_diss(micro, ls="none", marker="x", ax=plt.gca())

(argo_fine_depth.χ.sel(criteria="whalen") / 2).cf.plot.step(
    marker="o",
    color="C1",
    ax=plt.gca(),
    y="pressure_bins",
    label="$⟨χ⟩^{argo}_{Rω=3}/2$",
)

# (natre_fine_avg.χ.sel(criteria="whalen") / 2).cf.plot.step(
#    marker="o", color="C2", ax=plt.gca(), label="$⟨χ⟩_{fine-3}^{natre}/2$"
# )
(natre_fine_avg.χ.sel(criteria="whalen_7") / 2).cf.plot.step(
    marker="o", color="C3", ax=plt.gca(), label="$⟨χ⟩_{Rω=7}^{natre}/2$"
)

# (micro.chi / 2 - micro.Krho * micro.dTdz ** 2).cf.plot.step(
#    y="pres",
#   xscale="log",
#   xlim=(1e-11, 1e-8),
#   ylim=(1800, 200),
#   color="darkblue",
#   lw=1,
#   label="$⟨χ⟩/2 - K_ρ θ^m_z$",
# )

plt.legend(loc="right", bbox_to_anchor=(1.6, 0.5))
plt.xlabel("Variance production or dissipation [°C²/s]")
plt.gcf().set_size_inches((4.5, 6))
plt.title("NATRE")
plt.gcf().savefig("../images/natre-estimate.png")
_images/74582e3db4b6fcd676a9668315f037ec862be953c40c2dda4f759e62b21f9d7b.png

POP2 1/10° mesoscale variance budget#

Yiming is working with the following variance budget equation for the mesoscale after decomposing \(T = \bar{T} + T' + T'' ≡ T_m + T_e + T_t\) (the latter is Ferrari & Polzin notation). The \(T = \tilde{T} + T''\) decomposition is inherent when using cell-averaged equations.

(9)#\[\begin{equation} ∂_t \overline{T'^2} = -\overline{u'T'} ⋅ ∇_h\overline{T} - \overline{w'T'}∂_z\overline{T} - \overline{\vec{v} ⋅ ∇ T'^2/2} + \overline{T'∂_z \left[K ∂_zT'\right]} + … \end{equation}\]

Ferrari & Polzin’s equations (4.6, 4.7) for the mesoscale and microscale are (\(u\) here is the full 3D \(u\)):

(10)#\[\begin{align} ⟨u_eT_e⟩ ⋅ ∇T_m - ⟨\widetilde{u_t T_t} ⋅ ∇T_e⟩ &= 0 \ ⟨u_tT_t⟩ ⋅ ∇ T_m + ⟨\widetilde{u_t T_t} ⋅ ∇ T_e⟩ &= - \frac 12 ⟨χ⟩ \end{align}\]

Note that these equations are specialized for the “interior”.

For the model the averaging scales are:

  1. ~: grid-scale in vertical; single time-step in time i.e. the “microscale” is the “grid-scale”

  2. ⟨⟩: monthly in time; NATRE domain size in (x,y)

In Yiming’s notation these equations are

(11)#\[\begin{align} \overline{u'T'}⋅∇\overline{T} + \overline{w'T'}∂_z\overline{T} - \overline{T'∂_z \left[K ∂_zT' \right]} - \overline{T'A_h ∇^4 T'} = 0 \ K \left(∂_z(\overline{T} + T')\right)^2 &= - \frac 12 ⟨χ⟩ \end{align}\]

Now we can identify Yiming’s dissipation term (“VMIX” + “HDIFF”, ignoring solar heating etc.) as the scale transformation term in FP2005 and then use the microscale budget from FP2005 to write

(12)#\[\begin{equation} ⟨\widetilde{u_t T_t} ⋅ ∇T_e⟩ = - \frac 12 ⟨χ⟩ - ⟨u_tT_t⟩ ⋅ ∇ T_m \end{equation}\]

i.e. Yiming’s “mesoscale dissipation” term is really the difference between

  1. total microscale dissipation, and

  2. microscale production by stirring of the mean vertical gradient,

both of which can be estimated from microstructure observations.

Hide code cell source 
Image("../images/Tvariance_budget_NATREregion_with_EKE.png", width=600)
_images/596540b48690d1bbacd2fd65208fae64de30850faeec14e9a02c67eddbc29a27.png 
distributed.utils_perf - WARNING - full garbage collections took 13% CPU time recently (threshold: 10%)

Here are Yiming’s equations

Image("../yiming-budget.png")
_images/55429e08216b0b0404890cbdc07b41ea75e7c803841879823d3d138e6f899c38.png

POP 1/10 vs POP 1 vs ECCO#

TODO:#

  • [ ] Add ECCO IC

  • [ ] interpolate diffuivities to density surfaces?

  • [ ] Check pop-hires-annual-climatology.nc

  • [ ] For 1/10°, Infer 1° K_e, then average…?

  • [ ] Add Scott’s K; diagonal components; maybe compare cross-stream diffusivities?

  • [ ] Some tuning of Kappa for SMYLE.

  • [ ] 1/10° BGC run from MATT

  • [x] Email Keith about initialization for 1°.

  • [x] Compare ∇T for POP 1° and 1/10° initial conditions.

This plot compare POP simulations.

  • POP 1/10° is regridded to 1° grid.

  • All model output is averaged over years 42-61

  • 1° grid products are regridded to density space; then averaged in NATRE box for mean \(K_e, χ^{meso}\)

  • For \(|∇_ρT|\) I take 1° T field; and do plane fit in NATRE box; then estimate \(|∇T|^2\)

Conclusions:#

  1. POP 1/10° |∇T| agrees decently well with Argo |∇T| ; Argo and MIMOC agree quite well

  2. POP 1° initial conditions don’t capture the mean spice gradient in the NATRE region.

    • 1/10° initialized with quarter degree WOA?

    • 1° initialized with something else?

  3. POP 1° \(K\_{redi}\) is too high, below a 1000m, by a factor of 2-3.

  4. Inferring \(K_e = χ/|∇T|^2\) for the 1° gives basically the same value as KAPPA_ISOP

  5. Purely diffusive view of the problem. What if the production - med outflow - is too weak?

argo_delT2_plane.pres.attrs["positive"] = "down"
Hide code cell source 
plt.rcParams["figure.dpi"] = 140
plt.rcParams["figure.facecolor"] = "none"  # (0.91, 0.91, 0.91)
plt.rcParams["axes.facecolor"] = (0.99,) * 3
plt.rcParams.update(
    {
        "axes.spines.top": True,
        "axes.spines.bottom": True,
        "axes.spines.left": True,
        "axes.spines.right": True,
        "axes.edgecolor": (0.31,) * 3,
        "axes.grid.which": "both",
        "grid.linewidth": 0.5,
        "grid.alpha": 0.3,
        "grid.color": (0,) * 3,
    }
)

hirescolor = "darkorchid"
spinupcolor = "r"
spinupendcolor = "b"
lastcyclecolor = "deepskyblue"
color_ecco = "darkorange"

f, axx = plt.subplots(1, 3, sharey=True, constrained_layout=True)
ax = dict(zip(["K", "delT2", "chi"], axx))

legend_handles = []

for pop_, color, label in zip(
    [
        pop_natre["hires/clim"],
        # pop_1deg_spinup_profile.isel(decade=0),
        pop_natre["1deg/month1"],
        pop_natre["1deg/decadal"].isel(decade=5),
        # pop_1deg_natre_profile.sel(cycle=5),
    ],
    [hirescolor, spinupcolor, spinupendcolor, lastcyclecolor],
    ["POP2 1/10°", "POP2 1° month 1", "POP2 1° decade 6", "POP2 1° last cycle"],
):
    kwargs = dict(
        edges=ed.pop.get_edges(pop_.ds), orientation="horizontal", lw=2, color=color
    )
    ax["delT2"].stairs(np.sqrt(pop_["delT2_plane"].data), **kwargs)
    if "KAPPA_ISOP" in pop_:
        legend_handles.append(ax["K"].stairs(pop_["KAPPA_ISOP"].to_numpy(), **kwargs))
        ax["chi"].stairs(pop_["RediVar"].to_numpy(), **kwargs, label=label)

# K
legend_handles.append(
    groeskamp.Ke.cf.plot(
        y="pres", ax=ax["K"], color="limegreen", ls="--", label="Groeskamp (2020)"
    )
)
legend_handles.append(
    cole_natre.diffusivity.cf.plot(
        y="pres", ax=ax["K"], color="k", ls="--", label="Cole (2015)"
    )
)
legend_handles.append(ecco.Ke.cf.plot(y="vertical", ax=ax["K"], color=color_ecco))
# pop_hires_clim_profile.Ke.drop_vars("z_σ").cf.plot.step(
#    y="depth",
#    ax=ax["K"],
#    xlim=[0, 3000],
#    color=hirescolor,
#    label="POP 1/10° χ$^{e}$/|∇T|$^2$",
# )

# ∇T2
for data, color, label in zip(
    [
        argo_delT2_plane,  # mimoc_natre_profile.delT2_plane.rename({"pressure": "pres"})
        ecco_natre_profile.delT2_plane,
    ],
    [
        "k",
        color_ecco,  # "darkcyan"
    ],
    [
        "Argo",
        "ECCO Adjusted",  # "MIMOC"
    ],
):
    if label == "MIMOC":
        continue
    np.sqrt(data).cf.plot(
        y="vertical", ax=ax["delT2"], color=color, label=label, _labels=False
    )

# χmeso
(-1 * pop_natre["hires/clim"]["DISS"].drop_vars(["σ", "z_σ"])).cf.plot(
    ax=ax["chi"], lw=2, color=hirescolor, label="POP2 1/10°", _labels=False
)
(ecco.GM_CHI).sel(node="adj").cf.plot(
    ax=ax["chi"],
    y="vertical",
    lw=2,
    color=color_ecco,
    label=None,
    _labels=False,
)
# ecco.GM_CHI.sel(node="const").cf.plot(
#    ax=ax["chi"], y="vertical", lw=2, label="ECCO const", _labels=False
# )


dcpy.plots.errorbar(
    micro,
    x="residual",
    y="pres",
    color=mesomicrocolor,
    label="NATRE: $⟨χ⟩/2 - K_ρ^m ∂_zθ_m^2$",
    # lw=2.5,
    ax=ax["chi"],
    **error_kwargs,
)

ax["K"].set_title("")
ax["K"].set_ylim([3000, 200])
ax["K"].set_xlim([0, 2000])
ax["delT2"].set_xlim([1e-8, 4e-6])
ax["chi"].set_xlim([1e-12, 1e-8])
dcpy.plots.clean_axes(axx)

for a in axx:
    # a.legend(loc="upper center", bbox_to_anchor=(0.5, -0.0), frameon=False)
    a.xaxis.tick_top()
    a.xaxis.set_label_position("top")
    patch = mpl.patches.Rectangle(
        (0, 800), 1, 700, transform=a.get_yaxis_transform(), alpha=0.1, zorder=1
    )
    a.add_patch(patch)

ax["K"].set_xlabel("(a) $K_{e}$ [m$^2$/s]", labelpad=12)
ax["delT2"].set_xlabel("(b) $|∇_ρ^hT|^2$ [C/m]", labelpad=12)
ax["chi"].set_xlabel("(c) $χ_{e} = K_{e} |∇_ρ^hT|^2$ [C$^2$/s]", labelpad=12)

ax["K"].set_ylabel("depth [m]")
ax["chi"].set_xscale("log")
ax["chi"].xaxis.grid(True, which="both")
f.legend(loc="outside right")
# dcpy.plots.label_subplots(list(ax.values()), x=0.9, y=0.1)

f.set_size_inches((9, 6))
f.savefig("../images/coarse-models.pdf")
f.savefig("../images/coarse-models.png")
_images/5d1b9f10eb4aa693a1b8f20bb107bded1925b163be5a93cd0efb8a7ee33db645.png

ECCO#

  • Kux, Kvy are constant in time as expected.

  • But when I turn off the adjusted file, I still get values ≠ 1000

  • Interestingly adjusted \(K\) is higher right where we expect variance dissipation.

  • Adjust \(K\) is also high right where there is no isopycnal T gradient, so maybe overinterpreting?

TODO:

  • [x] compute isopycnal gradients

  • [ ] Make sure I understand the negative sign on Kuz, Kvz

  • [ ] compare χ from Cole to ECCO K adjusted map

  • [ ] I have trouble matching the published, remapped Kredi to Kux; issues might be the region I picked. Should get the data on the native grid.

kredi.sel(k=1000, method="nearest").sel(lat=slice(20, 40), lon=slice(290, 360)).hvplot()

(ecco.GM_Kux).cf.plot.line(hue="node", y="Z")
(ecco.GM_Kvy).cf.plot.line(hue="node", y="Z")
ecco.Ke.cf.plot(y="k", color="k")

[<matplotlib.lines.Line2D at 0x1783ec820>]
_images/a0047e3af7e27166381bca149f2d94d2e184a937648be191f494feacc10efbe9.png 
ecco.GM_CHI.cf.plot(hue="node", y="vertical")

[<matplotlib.lines.Line2D at 0x177eb07c0>,
 <matplotlib.lines.Line2D at 0x177f43a60>]
_images/62ece1a2ba911fc0cbb2580552570b77c8f19a586d9822d2ca23b2fe29caa7ea.png

3-term balance?#

pop_var = xr.open_zarr("../datasets/pop-variance-budget.zarr/")

pop_var

<xarray.Dataset>
Dimensions:  (depth: 62, lat: 2400, lon: 3600)
Coordinates:
    TLAT     (lat, lon) float64 dask.array<chunksize=(400, 800), meta=np.ndarray>
    TLONG    (lat, lon) float64 dask.array<chunksize=(400, 800), meta=np.ndarray>
  * depth    (depth) float64 0.05 0.15 0.25 0.35 ... 51.25 53.75 56.25 58.75
Dimensions without coordinates: lat, lon
Data variables:
    BC       (depth, lat, lon) float64 dask.array<chunksize=(62, 400, 800), meta=np.ndarray>
    HDIFF    (depth, lat, lon) float64 dask.array<chunksize=(62, 400, 800), meta=np.ndarray>
    PKC      (depth, lat, lon) float64 dask.array<chunksize=(62, 400, 800), meta=np.ndarray>
    VMIX     (depth, lat, lon) float64 dask.array<chunksize=(62, 400, 800), meta=np.ndarray>

SSH urms / Salinity on isopycnals#

argoclim = xr.open_zarr(
    "../datasets/argo-cole-clim-sigma2-2023-04-21.zarr/"
).cf.guess_coord_axis()
argoclim

<xarray.Dataset>
Dimensions:  (lat: 145, lon: 360, σ: 60)
Coordinates:
  * lat      (lat) float64 -64.5 -63.5 -62.5 -61.5 -60.5 ... 76.5 77.5 78.5 79.5
  * lon      (lon) float64 20.5 21.5 22.5 23.5 24.5 ... 376.5 377.5 378.5 379.5
  * σ        (σ) float32 34.15 34.15 34.17 34.18 34.22 ... 37.2 37.2 37.2 37.2
Data variables:
    CT       (lat, lon, σ) float64 dask.array<chunksize=(50, 50, 60), meta=np.ndarray>
    SA       (lat, lon, σ) float64 dask.array<chunksize=(50, 50, 60), meta=np.ndarray>
    z_σ      (lat, lon, σ) float32 dask.array<chunksize=(50, 50, 60), meta=np.ndarray>
hires = xr.open_zarr(
    "../datasets/pop-hires-annual-climatology-grid-vars.zarr",
    chunks={"nlat_t": 600, "nlon_t": 600},
)
hires["TLONG"] = xr.where(hires.TLONG > 180, hires.TLONG - 360, hires.TLONG)

/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/xarray/core/dataset.py:255: UserWarning: The specified Dask chunks separate the stored chunks along dimension "nlon_t" starting at index 600. This could degrade performance. Instead, consider rechunking after loading.
  warnings.warn(

import pop_tools

lores = xr.load_dataset("../datasets/pop-1deg-sigma2-42-61-clim.nc")

grid, logrid = pop_tools.to_xgcm_grid_dataset(pop_tools.get_grid("POP_gx1v7"))

lores["TLAT"] = grid.interp(logrid.ULAT, ["X", "Y"])
lores["TLONG"] = grid.interp(logrid.ULONG, ["X", "Y"])

lores["TLONG"] = xr.where(lores.TLONG > 180, lores.TLONG - 360, lores.TLONG)

OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import colorcet

kwargs = dict(
    y="latitude",
    x="longitude",
    vmin=35,
    vmax=36,
    add_colorbar=False,
)

f, ax = plt.subplots(
    1,
    4,
    constrained_layout=True,
    subplot_kw=dict(projection=ccrs.Orthographic(-20, 20)),
)


def plot_salt(da, ax):
    h = da.isel(σ=36).cf.plot(ax=ax, cmap="cet_bmy_r", **kwargs)
    da.isel(σ=36).cf.plot.contour(
        ax=ax, colors="darkgray", linewidths=1, levels=np.linspace(35, 36, 11), **kwargs
    )
    return h


plot_salt(argoclim.SA.sel(lat=slice(10, 50), lon=slice(360 - 40, 360)), ax[0])
plot_salt(hires.SALT.isel(nlat_t=slice(1250, 1650), nlon_t=slice(700, 1140)), ax[1])
plot_salt(lores.SALT.cf.isel(X=slice(0, 35), Y=slice(225, 320)), ax[2])
h = plot_salt(ecco_tile2.SALT, ax=ax[3])

[ed.plot.make_nice_natre_map(ax_) for ax_ in ax]

ax[0].set_title("(a) Argo")
ax[0].gridlines(draw_labels=["left", "bottom"])

ax[1].set_title("(b) POP 1/10°")
ax[1].gridlines(draw_labels=["bottom"])

ax[2].set_title("(c) POP 1°")
ax[2].gridlines(draw_labels=["bottom"])

ax[3].set_title("(d) ECCO")
ax[3].gridlines(draw_labels=["right", "bottom"])

f.colorbar(
    h,
    ax=ax,
    orientation="horizontal",
    extend="both",
    label="Salinity on $σ_{2}$=36.33 kg/m3",
    aspect=40,
    shrink=0.6,
)

f.set_size_inches((10, 4))

f.savefig("../images/isopycnal-salt.png", dpi=200, bbox_inches="tight")
_images/2927dbaf6895a11a689b886716e629df8bb1d2c0ebe7d1b2e4ac6d16f8ac30a2.png 
2023-04-26 22:15:52,337 - tornado.application - ERROR - Exception in callback <bound method SystemMonitor.update of <SystemMonitor: cpu: 1 memory: 15 MB fds: 33>>
Traceback (most recent call last):
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/tornado/ioloop.py", line 921, in _run
    val = self.callback()
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/system_monitor.py", line 161, in update
    net_ioc = psutil.net_io_counters()
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/psutil/__init__.py", line 2114, in net_io_counters
    rawdict = _psplatform.net_io_counters()
OSError: [Errno 12] Cannot allocate memory
2023-04-27 08:50:02,100 - distributed.scheduler - WARNING - Worker failed to heartbeat within 300 seconds. Closing: <WorkerState 'tcp://127.0.0.1:59063', name: 2, status: running, memory: 0, processing: 0>
2023-04-27 08:50:02,330 - distributed.scheduler - WARNING - Worker failed to heartbeat within 300 seconds. Closing: <WorkerState 'tcp://127.0.0.1:59065', name: 0, status: running, memory: 0, processing: 0>
2023-04-27 08:50:02,417 - distributed.scheduler - WARNING - Received heartbeat from unregistered worker 'tcp://127.0.0.1:59065'.
2023-04-27 08:50:02,429 - distributed.scheduler - WARNING - Received heartbeat from unregistered worker 'tcp://127.0.0.1:59063'.
2023-04-27 08:50:02,436 - distributed.worker - ERROR - Scheduler was unaware of this worker 'tcp://127.0.0.1:59063'. Shutting down.
2023-04-27 08:50:02,434 - distributed.worker - ERROR - Scheduler was unaware of this worker 'tcp://127.0.0.1:59065'. Shutting down.

Traceback (most recent call last):
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/utils.py", line 752, in wrapper
    return await func(*args, **kwargs)
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/worker.py", line 1525, in close
    await self.finished()
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/core.py", line 502, in finished
    await self._event_finished.wait()
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/asyncio/locks.py", line 226, in wait
    await fut
asyncio.exceptions.CancelledError

Traceback (most recent call last):
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/utils.py", line 752, in wrapper
    return await func(*args, **kwargs)
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/worker.py", line 1525, in close
    await self.finished()
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/core.py", line 502, in finished
    await self._event_finished.wait()
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/asyncio/locks.py", line 226, in wait
    await fut
asyncio.exceptions.CancelledError
2023-04-27 08:50:06,666 - distributed.nanny - ERROR - Worker process died unexpectedly
Exception in thread Nanny stop queue watch:
Traceback (most recent call last):
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/threading.py", line 973, in _bootstrap_inner
2023-04-27 08:50:06,679 - distributed.nanny - ERROR - Worker process died unexpectedly
2023-04-27 08:50:08,541 - distributed.nanny - WARNING - Restarting worker
2023-04-27 08:50:08,693 - distributed.nanny - WARNING - Restarting worker
2023-04-27 08:50:29,898 - distributed.nanny - WARNING - Worker process still alive after 3.1999897766113286 seconds, killing
2023-04-27 08:50:29,963 - distributed.nanny - WARNING - Worker process still alive after 3.199994201660157 seconds, killing
2023-04-27 08:50:29,998 - distributed.nanny - ERROR - Error in Nanny killing Worker subprocess
Traceback (most recent call last):
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/asyncio/tasks.py", line 452, in wait_for
    fut.result()
asyncio.exceptions.CancelledError

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/nanny.py", line 607, in close
    await self.kill(timeout=timeout, reason=reason)
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/nanny.py", line 390, in kill
    await self.process.kill(reason=reason, timeout=0.8 * (deadline - time()))
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/nanny.py", line 844, in kill
    await process.join(max(0, deadline - time()))
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/process.py", line 330, in join
    await wait_for(asyncio.shield(self._exit_future), timeout)
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/utils.py", line 1849, in wait_for
    return await asyncio.wait_for(fut, timeout)
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/asyncio/tasks.py", line 454, in wait_for
    raise exceptions.TimeoutError() from exc
asyncio.exceptions.TimeoutError
2023-04-27 08:50:30,100 - distributed.nanny - ERROR - Error in Nanny killing Worker subprocess
Traceback (most recent call last):
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/asyncio/tasks.py", line 452, in wait_for
    fut.result()
asyncio.exceptions.CancelledError

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/nanny.py", line 607, in close
    await self.kill(timeout=timeout, reason=reason)
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/nanny.py", line 390, in kill
    await self.process.kill(reason=reason, timeout=0.8 * (deadline - time()))
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/nanny.py", line 844, in kill
    await process.join(max(0, deadline - time()))
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/process.py", line 330, in join
    await wait_for(asyncio.shield(self._exit_future), timeout)
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/distributed/utils.py", line 1849, in wait_for
    return await asyncio.wait_for(fut, timeout)
  File "/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/asyncio/tasks.py", line 454, in wait_for
    raise exceptions.TimeoutError() from exc
asyncio.exceptions.TimeoutError

aviso = xr.open_zarr("~/datasets/aviso_monthly.zarr", consolidated=True)
aviso["longitude"] = aviso.longitude - 360
aviso = aviso.sel(latitude=slice(10, 50), longitude=slice(-40, 5))
aviso["urms"] = np.sqrt(aviso.ugosa**2 + aviso.vgosa**2)
aviso["urms_mean"] = aviso.urms.mean("time").load()
aviso.urms_mean

<xarray.Dataset>
Dimensions:    (time: 312, latitude: 160, longitude: 160, nv: 2)
Coordinates:
  * latitude   (latitude) float32 10.12 10.38 10.62 10.88 ... 49.38 49.62 49.88
  * longitude  (longitude) float32 -39.88 -39.62 -39.38 ... -0.625 -0.375 -0.125
  * nv         (nv) int32 0 1
  * time       (time) datetime64[ns] 1993-01-01 1993-02-01 ... 2018-12-16
Data variables:
    adt        (time, latitude, longitude) float64 dask.array<chunksize=(12, 160, 160), meta=np.ndarray>
    crs        int32 ...
    lat_bnds   (latitude, nv) float32 dask.array<chunksize=(160, 2), meta=np.ndarray>
    lon_bnds   (longitude, nv) float32 dask.array<chunksize=(160, 2), meta=np.ndarray>
    sla        (time, latitude, longitude) float64 dask.array<chunksize=(12, 160, 160), meta=np.ndarray>
    ugos       (time, latitude, longitude) float64 dask.array<chunksize=(12, 160, 160), meta=np.ndarray>
    ugosa      (time, latitude, longitude) float64 dask.array<chunksize=(12, 160, 160), meta=np.ndarray>
    vgos       (time, latitude, longitude) float64 dask.array<chunksize=(12, 160, 160), meta=np.ndarray>
    vgosa      (time, latitude, longitude) float64 dask.array<chunksize=(12, 160, 160), meta=np.ndarray>
    urms       (time, latitude, longitude) float64 dask.array<chunksize=(12, 160, 160), meta=np.ndarray>
    urms_mean  (latitude, longitude) float64 0.1002 0.1002 ... 0.1431 0.1444
Attributes: (12/45)
    Conventions:                     CF-1.6
    Metadata_Conventions:            Unidata Dataset Discovery v1.0
    NCO:                             netCDF Operators version 4.7.9 (Homepage...
    cdm_data_type:                   Grid
    comment:                         Sea Surface Height measured by Altimetry...
    contact:                         servicedesk.cmems@mercator-ocean.eu
    ...                              ...
    summary:                         SSALTO/DUACS Delayed-Time Level-4 sea su...
    time_coverage_duration:          P1D
    time_coverage_end:               1993-01-01T12:00:00Z
    time_coverage_resolution:        P1D
    time_coverage_start:             1992-12-31T12:00:00Z
    title:                           DT merged all satellites Global Ocean Gr...
import cartopy.crs as ccrs

f, ax = plt.subplots(
    1,
    2,
    constrained_layout=True,
    subplot_kw=dict(projection=ccrs.Orthographic(-20, 20)),
)
aviso.urms_mean.plot(
    transform=ccrs.PlateCarree(),
    robust=True,
    cmap="Reds",
    vmin=0,
    ax=ax[0],
    add_colorbar=False,
    vmax=0.15,
)
ed.plot.make_nice_natre_map(ax[0])
ax[0].set_title("(a) AVISO")
ax[0].gridlines(draw_labels=["left", "bottom"])

<cartopy.mpl.gridliner.Gridliner at 0x16df1a790>
_images/a313bb37477983db92f51c5e3a13dacf44fb5332d7281078dcd19f5f6f1d9980.png

Vertical grid spacings#

ECCO#

ecco.drC.set_xindex("Zp1").sel(Zp1=slice(800, 1500))

<xarray.DataArray 'drC' (k_p1: 7)>
array([ 95.27 ,  97.415,  98.75 ,  99.63 , 100.67 , 102.945, 107.945],
      dtype=float32)
Coordinates:
  * k_p1     (k_p1) int64 27 28 29 30 31 32 33
    face     int64 2
  * Zp1      (k_p1) float32 861.5 958.0 1.056e+03 ... 1.357e+03 1.461e+03
    drC      (k_p1) float32 95.27 97.42 98.75 99.63 100.7 102.9 107.9
    PHrefF   (k_p1) float32 8.451e+03 9.398e+03 ... 1.331e+04 1.434e+04
Attributes:
    standard_name:  cell_z_size_at_w_location
    long_name:      cell z size
    units:          m

POP 1/10°#

pop_natre["hires/clim"].sel(depth=slice(800, 1500))["depth"].diff("depth")

<xarray.DataArray 'depth' (depth: 4)>
array([105.88061523, 121.49841309, 138.3626709 , 155.93029785])
Coordinates:
  * depth    (depth) float64 984.7 1.106e+03 1.245e+03 1.4e+03
    time     object 0052-01-01 05:17:48.750000
    σ        (depth) float64 36.57 36.68 36.77 36.86
    z_σ      (depth) float64 984.7 1.106e+03 1.245e+03 1.4e+03
Attributes:
    positive:       down
    standard_name:  depth
    units:          meter

POP 1°#

pop_natre["1deg"]

<xarray.DatasetView>
Dimensions:  ()
Data variables:
    *empty*

Error bounds#

Microstructure#

Ferrari & Polzin (2005):

image.png

Studying Cole estimate#

Previously I was seeing some big disagreement. This was due to gradients estimated in depth-space instead of in isopycnal space.

Compare Cole & Groeskamp diffusivities#

Reproducing Groeskamp Figure 1. This seems to work well; but note the different domain. I saw this in his MATLAB script.

%% MIXING IN NATRE REGION
% THE NATRE REGIONS LIES BETWEEN
% 319.5 - 344.5 degrees lon
%  19.5 -  29.5 degrees lat
NATRE_x = 320:345;
NATRE_y = 110:120;

cole = ed.read_cole(resolution="3deg").sel(
    lat=slice(19.5, 29.5), lon=slice(319.5, 344.5)
)
groeskamp = (
    xr.open_dataset("/home/deepak/datasets/groeskamp2020/groeskamp2020.nc")
    .sel(lat=slice(19.5, 29.5), lon=slice(319.5, 344.5))
    .rename({"depth": "pres"})
)


def plot_mean_std(da, label):
    mean = da.mean(["lat", "lon"])
    std = da.std(["lat", "lon"])
    (hdl,) = mean.cf.plot(ylim=(2000, 0), label=label)
    dcpy.plots.fill_between(
        xr.concat((mean - std, mean + std), dim="err"),
        axis="x",
        x="err",
        y="pres",
        alpha=0.3,
    )


plot_mean_std(groeskamp.Ke_0, label="Groeskamp unsuppressed")
plot_mean_std(groeskamp.Ke, label="Groeskamp suppressed")
plot_mean_std(cole.diffusivity, label="Cole")
plt.gcf().set_size_inches((3.5, 5))
plt.gca().set_xlim((-500, 3000))
plt.legend(loc="center left", bbox_to_anchor=(1, 0.5))

/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice.
  var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof,

<matplotlib.legend.Legend at 0x7f6a3f8c57f0>
_images/38fe3a99bcbdd55646b62077ee1199a362c75f8d958ad71482461c6584d8499c.png

Depth-vs-density space estimation of \(K_e |∇T|²\)#

Turns out this does not matter; I was worried that diffusivity in depth space and |∇T| estimated in density space might combine badly. But moving diffusivity to sigma space or |∇T| to depth space results in same mean variance tendency

cole = ed.read_cole()

import xgcm

cole["diffusivity_sig"] = xgcm.transform.linear_interpolation(
    cole.diffusivity, cole.density_mean_depth, cole.sigma, "pres", "pres", "sigma"
)

# read, take mean in region; swap to pressure coordinate
colegrad = xr.open_dataset("../datasets/cole-clim-gradient-natre.nc").sel(
    lat=slice(24, 28), lon=slice(360 - 30, 360 - 25)
)
colegrad["sigma"] = cole.sigma

cole_natre = cole[
    ["diffusivity", "diffusivity_sig", "density_mean_depth", "depth_mean_sig"]
].interp(lat=colegrad.lat.data, lon=colegrad.lon.data)
dTiso_pres = xgcm.transform.linear_interpolation(
    colegrad.dTiso,
    colegrad.pres,
    cole.pres,
    "sigma",
    "sigma",
    "pres",
)

cole_var_sig = (
    (cole_natre.diffusivity_sig * colegrad.dTiso**2)
    .mean(["lat", "lon"])
    .assign_coords(pres=cole_natre.depth_mean_sig.mean(["lat", "lon"]))
)
cole_var_pres = (cole_natre.diffusivity * dTiso_pres**2).mean(["lat", "lon"])

cole_var_sig.plot(y="pres")
cole_var_pres.plot(y="pres", yincrease=False)
plt.legend(["sigma space", "depth space"])

<matplotlib.legend.Legend at 0x7fa53eb855b0>
_images/54737be689897ce406c7fccaf3c803f414d9b4e50177261d9370c6dda1d9946e.png

Depth space calculation#

The following is when I was using \(|∇T|\) estimated in depth space and multiplying by a \(K_e\) provided in depth space

Initial surprising observation: We expect to see peak \(K_e\) between 600 and 1600m but we mostly don’t! Both \(⟨S'S'⟩\) and \(|∇S|\) are large in the depth we expect a peak. There seems to be approximate agreement in the 1000-1200m range at (-34°E, -32°E) but missing data at deeper depths

Here are estimates of variance produced by eddy stirring at points in the NATRE box from Cole et al (2015). There is some spread but the peak between 400 and 1000 seems robust. Thick black line is the mean over the box.

cole_var_all.stack(latlon=["lat", "lon"]).drop("latlon").cf.plot.line(
    y="pres", add_legend=False, size=4, aspect=0.8, xscale="log"
)
cole_var.cf.plot(color="k", lw=2)
plt.xlabel("$K_e^{cole} |∇T_{iso}^{argo}|²$ [°C²/s]")
# (cole.diffusivity.interp_like(dTiso) * colegrad.dSiso**2).mean(["lat", "lon"]).cf.plot()

Text(0.5, 0, '$K_e^{cole} |∇T_{iso}^{argo}|²$ [°C²/s]')
_images/5a2fafd95cb62c71c37708fa39821d062996780f3567e3621142857bb3640c99.png

These are the various terms in the Cole analysis. Note minimum in ∇S between 200 and 800, and associated high mixing lengths ≳ 300km.

ed.plot.plot_cole_profile(lat=26, lon=360 + np.array([-35, -34, -33, -32, -30, -27]))
_images/980f86cd5f2e76ed3639718b4fee30f3f699b82f6c8ab76fefbe335797ad8114.png

T-S diagrams: Argo vs. NATRE#

Lets compare T-S to see if they are similar. I am exploring ways to

  1. interpret the “high” variance production values between 400m and 800m for the Cole estimate. Maybe the NATRE T-S is not representative of the argo T-S?

    • It is representative!

  2. explain why the Argo ⟨χ⟩ seems to miss any eddy contribution. Is the scatter less at those density levels?

    • There is some undersampling (of course) because the sampling dz is low below 1000m usually.

    • It seems to capture some of the spread; but the NATRE dataset seems to fill up the space more.

The density contours are ρ corresponding to depths np.arange(200, 2000, 200) in the Cole dataset. Clearly the top few contours see not much T-S stirring and the T-S relationship is “tight”. Ferrari & Polzin interpreted this as a balance between microscale stirring and dissipation.

  • So the high K_e estimate between 400m and 800m in Cole et al (2015) is not supported by the TS spread.

  • Though interestingly the Groeskamp estimate has the same issue then.

Is the finestructure estimate undersampling eddy stirring periods? Compare T-S spread for

  1. all argo profiles

  2. argo profiles where I can make a finestructure estimate.

  3. NATRE profiles

Hide code cell source 
argo_natre = xr.load_dataset("../datasets/argo/natre.nc")
argo_natre["THETA"] = dcpy.eos.ptmp(
    argo_natre.PSAL, argo_natre.TEMP, argo_natre.PRES, natre.reference_pressure.item()
)
argo_fine_full_profiles = xr.load_dataset("argo_profiles_used_for_finestructure.nc")

# choose density levels that highlight the depth range we are interested in
rho_levels = (
    cole.density_mean_depth.interp_like(colegrad)
    .mean(["lat", "lon"])
    .sel(pres=np.arange(200, 2000, 200), method="nearest")
)

All argo vs NATRE#

Hide code cell source 
kwargs = dict(
    hexbin=False,
    size=0.5,
    plot_kwargs={"alpha": 0.1},
    Sbins=np.arange(35, 36, 0.05),
    Tbins=np.arange(2, 15, 0.5),
    Pref=0,
)

_, ax = dcpy.oceans.TSplot(
    argo_natre.PSAL,
    argo_natre.THETA,
    **kwargs,
    rho_levels=rho_levels,
)

# mask = argo_fine_full_profiles.fine_deepest_bin > 1000

# _, ax = dcpy.oceans.TSplot(
#    argo_fine_full_profiles.where(mask, drop=True).PSAL,
#    argo_fine_full_profiles.where(mask, drop=True).THETA,
#    **kwargs,
#    rho_levels=rho_levels,
# )

dcpy.oceans.TSplot(
    natre.salt,
    natre.theta,
    ax=ax,
    **kwargs,
)

ax["s"].legend(["argo full", "natre"])
plt.gcf().savefig("../images/ts-natre.png")
_images/012e42d5095eff0798b5ec595097c6804096fb9559a22463c595e3f630af73c9.png

Finestructure Argo vs NATRE#

Hide code cell source 
kwargs = dict(
    hexbin=False,
    size=0.5,
    plot_kwargs={"alpha": 0.1},
    Sbins=np.arange(35, 36, 0.05),
    Tbins=np.arange(2, 10, 0.5),
    Pref=0,
)

The distributions look similar and the argo data capture more extremes; but the NATRE data fill a lot of the space (though smaller S.D.). Maybe the average isn’t that great and it’s totally undersampling

Hide code cell source 
mask = argo_fine_full_profiles.fine_deepest_bin > 1000

_, ax = dcpy.oceans.TSplot(
    natre.salt,
    natre.theta,
    **kwargs,
    rho_levels=rho_levels,
)

dcpy.oceans.TSplot(
    argo_fine_full_profiles.where(mask, drop=True).PSAL,
    argo_fine_full_profiles.where(mask, drop=True).THETA,
    ax=ax,
    **kwargs,
)

ax["s"].legend(["natre", "argo fine"])
plt.gcf().savefig("../images/ts-natre.png")
_images/2d6a61598d8d6e2394221c02655c3332b467f6f876d8e4add3f84e0a8db030ef.png

It’s definitely a big undersampling when compared to all Argo profiles in the area. How good is good enough?

Hide code cell source 
_, ax = dcpy.oceans.TSplot(
    argo_natre.PSAL,
    argo_natre.THETA,
    **kwargs,
    rho_levels=rho_levels,
)

mask = argo_fine_full_profiles.fine_deepest_bin > 1000


masked = argo_fine_full_profiles.where(mask, drop=True)
_, ax = dcpy.oceans.TSplot(
    masked.PSAL,
    masked.THETA,
    ax=ax,
    rho_levels=None,
    **kwargs,
)


masked.groupby
ax["s"].legend(["argo full", "argo fine"])
plt.gcf().savefig("../images/ts-natre-argo.png")
_images/2971b8da803d7294c3ba9b6320bf714bc81b4b2472a4047d661fc46180596916.png

Conclusions#

  1. Can we define a metric using (T-S spread in all argo profiles) vs (TS spread in profiles where we can make a finestructure estimate) to decide where we might be able to make an estimate using Argo?

  2. Can we refine this?

    • If we expect regions above 1000m then it might not be so bad.

    • Could make a bigger box? (skeptical)

  3. This really points out the limitation of the finestructure estimate, at least deeper down?

    • How does the Whalen scatter plot work out?

Is there spatial uncertainty? Cole estimate: Are we looking in the wrong place? No#

Here’s the NATRE region marked on a map of \(K_e\). There is just low diffusivity everywhere in the box; but may be the mean signal is just northwest of the box

fg = (
    cole.diffusivity.sel(pres=[800, 1200], method="nearest")
    .sel(lon=slice(300, 350), lat=slice(20, 40))
    .plot(robust=True, row="pres", cbar_kwargs={"orientation": "horizontal"})
)

fg.map(
    lambda: plt.plot(
        360 + np.array([-31, -26, -26, -31, -31]), [24, 24, 28, 28, 24], lw=2
    )
)

<xarray.plot.facetgrid.FacetGrid at 0x7fbb1793f2b0>
_images/d259c44da201fda2581f3d4dfeadb781a8c652e9c51232e52edb21fbc94c1d5d.png